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Abstract: In order to make the prediction of land surface heat fluxes more robust, two 
improvements were made to an operational two-layer model proposed previously by Zhang. 
These improvements are: 1) a surface energy balance method is used to determine the 
theoretical boundary lines (namely 'true wet/cool edge' and 'true dry/warm edge' in the 
trapezoid) in the scatter plot for the surface temperature versus the fractional vegetation 
cover in mixed pixels; 2) a new assumption that the slope of the T m - f curves is mainly 
controlled by soil water content is introduced. The variables required by the improved 
method include near surface vapor pressure, air temperature, surface resistance, 
aerodynamic resistance, fractional vegetation cover, surface temperature and net radiation. 
The model predictions from the improved model were assessed in this study by in situ 
measurements, which show that the total latent heat flux from the soil and vegetation are in 
close agreement with the in situ measurement with an RMSE (Root Mean Square Error) 
ranging from 30 w/m ~50 w/m 2 , which is consistent with the site scale measurement of 
latent heat flux. Because soil evaporation and vegetation transpiration are not measured 
separately from the field site, in situ measured CO2 flux is used to examine the modeled 
/LEVeg. Similar trends of seasonal variations of vegetation were found for the canopy 
transpiration retrievals and in situ CO2 flux measurements. The above differences are 
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mainly caused by 1) the scale disparity between the field measurement and the MODIS 
observation; 2) the non-closure problem of the surface energy balance from the surface 
fluxes observations themselves. The improved method was successfully used to predict the 
component surface heat fluxes from the soil and vegetation and it provides a promising 
approach to study the canopy transpiration and the soil evaporation quantitatively during the 
rapid growing season of winter wheat in northern China. 

Keywords: Two-layer model, surface evapotranspiration, surface energy balance, Bowen 
Radio, trapezoid method 



1. Introduction 

Evapotranspiration (IE, soil evaporation and vegetation transpiration) from the land surface is an 
important link between the surface energy balance and the hydrologic cycle. Its accurate 
characterization is therefore very important in the study of the terrestrial ecosystem, climate dynamics 
and hydrologic cycle. At present, estimate of regional evapotranspiration has been made possible by 
using the remote sensing observations in combination with the surface meteorological data. In the past 
years, several remote sensing methods were developed to simulate surface-atmosphere interactions and 
to retrieve the terrestrial evapotranspiration over a wide range of spatial scales [1]. 

By treating the soil-vegetation system as a single uniform leaf, the big-leaf model simplified the 
mechanism of the energy exchange between the surface and the atmosphere, and therefore the regional 
scale evapotranspiration simulation is made. This category of models is simple and convenient to use, 
but the limitation is that this big-leaf approximation in the model is not applicable to surfaces with 
highly spatial heterogeneity due to large differences of surface energy exchange between soil and 
vegetation, such as in arid or semi-arid areas. Therefore, a two-layer model is proposed and the surface 
available energy is partitioned between soil and vegetation to overcome the limitation of the big-leaf 
model. These models have an improvement over the big-leaf models when applied to sparsely 
vegetated surfaces [2]. Another motivation to study the two-layer model is that China exhibits a highly 
heterogeneous land cover due to a large population and scattered residential areas. The two-layer 
model is more suitable in such areas. 

In the existing two-layer models, the cores of the algorithm primarily lie in two aspects: 
(1) accurately decomposing surface temperature of mixed pixel (T m ) into soil temperature (T so n) and 
vegetation temperature (T veg ); (2) obtaining accurate surface resistances, such as aerodynamic 
resistance, canopy resistance, residual resistance. In recent years, many attempts had been made to 
investigate the two issues. For instance, Norman and Kustas [3,4,5] used remote measurements of 
surface directional brightness temperature and some ancillary data to obtain soil temperature and 
vegetation temperature (called multi-angle method), applied Beer's law to partition net radiation of 
mixed pixel (R n ) and employed Monin-Obukhov similarity theory to compute aerodynamic resistance 
(r a ); Zhang et al. [6] presented a PCACA algorithm (PCACA, Pixel Component Arranging and 
Comparing Algorithm) and a layered energy-separating algorithm on the basis of triangle method and 
Bowen-ratio energy balance method to partition surface temperature, surface albedo (oc m ) and net 
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radiation of mixed pixel, and finally to estimate soil evaporation (AEsoii) and vegetation transpiration 
(/LEVeg). Because the multi-angle satellite data is not always available, multi-angle method of surface 
temperature decomposing is limited for applications, In contrast, PCACA algorithm is more 
convenient because only single angle remote sensed data are required and it can be provided from most 
of the satellite data. Additionally, by using the layered energy-separating algorithm the core of which 
is Bowen-ratio energy balance method, the uncertainties in surface energy partitioning based on the 
Beer's law are reduced. 

PCACA algorithm and layered energy-separating algorithm utilize the scatter plot of the surface 
temperature against vegetation fraction cover (T m -f space) to determine the soil water status, like the 
trapezoid method. On the basis of: (1) the assumption that the configuration of T m - f space is not 
primarily caused by differences in atmospheric conditions and soil attributes (eg. air temperature T a , 
aerodynamic resistance r a , surface reflectivity a) but by the variations of soil water availability; (2) the 
fact that iso-lines of equal soil water availability are nearly straight in T m -f space, which was reported 
in previous studies on trapezoid method [7-11], Zhang et al. [6] indicated that T so i\ values for all pixels 
at an iso-line are equivalent, so are for T veg values, which is just like the case that while measuring the 
same area constituted by soil and vegetation at varying view angles, T so n and r veg are invariable and 
thus T m observations only vary with vegetation fraction cover. Under this condition, T so n and r veg could 
be obtained by calculating the slopes of these iso-lines of equal soil water availability (dTJdf). 
Detailed descriptions about the two algorithms and the trapezoid method can be found in Zhang et al. [6] 
and in Carlson et al. [9, 12], respectively. In this approach, an assumption of a uniform atmospheric 
environment and homogeneous soil surface is required. In most cases, however, this strict requirement 
can't be satisfied, especially on a regional scale. In addition, the identification of the trapezoid shape of 
T m -f space in the trapezoid method requires a flat surface and a large number of pixels over an area 
with a wide range of soil wetness and vegetation fraction cover, which usually cannot be satisfied 
within a limited study area, and therefore some subjectivity would be introduced in the determination 
of the trapezoid shape and it will inevitably introduce some extra errors in the T soi \ and T veg 
calculations, The temperature separation finally influences XE estimate. 

In this paper, to improve the accuracy ofXE estimate using the two-layer model presented by Zhang 
et al. [6], two modifications were made to the PCACA algorithm and layered energy-separating 
algorithm mentioned above: (1) to ensure that the assumption that the configuration of T m -f space is 
mainly controlled by soil water availability is reasonable, the effects of atmospheric conditions on 
surface temperature would be ignored by using the averaging method; (2) to identify the shape of the 
trapezoid bounded by true wet/cool edge' and 'true dry/warm edge, a general method based on surface 
energy balance was used. Finally, comparisons between XE retrievals from the original model and the 
improved one, and in situ XE measurements were done to assess the improved method. 

2. Model Descriptions 

The two-layer model used in this study was presented by Zhang et al. [6], and in it the PCACA 
algorithm and layered energy-separating algorithm are the key algorithms. In the model, land cover is 
simplified as a mixture of two elements, namely, vegetation and bare soil. The energy fluxes are 
partitioned between the soil and vegetation, and energy exchange between vegetation and bare soil is 
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negligible. Two parameters, albedo and surface temperature, are the main different characteristics of 
vegetation and bare soil, and lead to different interactions between them and atmosphere. Air 
temperature, air humidity, wind speed and aerodynamic resistance are approximated the same for 
vegetation and bare soil in the same pixel due to intensively atmospheric blending effect. 

2. 1 An interpretation of T m -f space 

Figure 1 provides a conceptual illustration of T m - f space, where "true wet/cool edge" of the 
trapezoid is related to surface conditions of potential evapotranspiration and has minimum surface 
resistance to evapotranspiration (r smin ). 

Figure 1. Scatter plot of surface temperature against vegetation fraction cover true dry edge 
observed dry edge observed wet edge true wet edge iso-line of equal vegetation fraction cover 
iso-line of equal soil water availability. 




0.0 0.30 0.60 0.90 1.20 

Fractional vegetation cover 



Soil water content equals to field capacity; 'true dry/warm edge' represents zero evapotranspiration 
and has maximum surface resistance to evapotranspiration (r smax ). If the positions of the two edges are 
determined, the shape and the structure of the trapezoid can be fixed and the consequent calculations 
of the surface heat fluxes could be done. To make the illustration easier to follow, four points are 
defined: r s d and r sw represent the points of true dry bare soil and true water saturated bare soil, 
respectively. 7Vd represents true dry full-cover vegetation and T vw represents true water saturated, full- 
cover vegetation. The above definitions all correspond to the ideal surface conditions, namely there 
exist driest and wettest bare soil and full-cover vegetation. In reality there are always insufficient 
number of pixels that can cover all kinds of soil wetness and vegetation fraction cover within the study 
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area, which leads to a difficulty in determining "true wet/cool edge" and "true dry/warm edge", as a 
result, "observed wet/cool edge" and "observed dry/warm edge" (dashed lines) are often defined 
according to the envelop shape of the actual scatter plot to represent actual two extreme soil moisture 
conditions and are used to replace "true wet/cool edge" and "true dry/warm edge" in most applications, 
although some errors would be introduced. Iso-line of equal vegetation fraction cover intersects "true 
dry edge" and "true wet edge" at true maximum temperature and true minimum temperature denoted 
as r m i jmax and T"mi,min, respectively. It has to be noted that for the trapezoid constructed by data with 
coarser pixel resolution (eg. 1km), surface temperature at "true dry edge" generally is higher than that 
at "observed dry edge", contrarily surface temperature at "true wet edge" is lower than that at 
"observed wet edge" according to the findings of Carlson [9] about the scale issues of trapezoid space. 
The major reason is the absence of the driest and wettest pixels due to the great mixing effect of the 
pixel (vegetation/bare soil). 

2.2. PCACA Algorithm 

The PCACA algorithm is a method to decompose the mixed surface temperature. As mentioned 
above, its physical basis is the fact that soil water status can be represented by the configuration of T m - 
f space, and like the triangle method, the crucial assumption is that surface temperature is mainly 
controlled by soil water availability [9, 10, 12]. 

The basic formulations used to decompose mixed surface temperature are shown in Eq. (1) and Eq. 
(2). The derivations of Eq. (2) can be found in Appendix I. 

^m T m= (T£ V egf r ^g +(T£ soil( l -f) T Ll 0) 

k = ~ T veg - T soil (2) 

where, e m , e veg and e so ii are the broadband emissivities of mixed pixel, vegetation and bare soil, 
respectively; a is the Stefan-Boltzmann constant; dr m /d/ represents the slope of iso-line of equal soil 
water availability, expressed as k in the remainder of the paper. In the study, constant e veg of 0.97 and 
£ so ii of 0.95 were used. By simply weighting the fractional cover for vegetation and bare soil, mixed 
pixel emissivity e m for each pixel can be computed following Eq. (3) [13]. 

£ m =f £ ve g + Q-f) £ soil ( 3 ) 

From the above three equations, we can see that to solve T &0 [\ and r veg , k for each pixel is the key 
parameter and needs to be obtained firstly. 

According to Figure 1, two procedures were needed to derive k in the study: (1) determining the 
shape of the trapezoid bounded by 'true dry edge' and 'true wet edge', namely locating the 'true dry 
edge' and 'true wet edge' on the trapezoid, and thereby calculating the two slopes of the two bounder 
lines (k u is for upper bound, At is for lower bound). A detailed method for this will be illustrated in 
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section 3; (2) linearly interpolating k between the highest and the lowest surface temperature to obtain 
the slope for each pixel based on the equation below: 

T -T 



mi mi.mra 
m;, max mi,mm 



(k u -k L ) + k L (4) 



Considering that the precise formulation of the relationship is actually unknown, linear 
interpolation is a reasonable approximation according to the study on the configuration of T m -f space 
by Calson [9, 10] and Moran [14]. 

2.3 Layered Energy-separating Algorithm 

Essentially, the aim of the Layered Energy-separating algorithm is to calculate Bowen-ratio (fi, the 
radio of sensible heat flux to latent heat flux) of soil and vegetation expressed as /? so ii and /? veg , 
respectively. By using the relationship between Water Deficit Index (WDT), evapotranspiration and 
potential evapotranspiration (XEq) illustrated by Moran et al. [14], Eq. (5) can be derived: 

T -T 

P ^ m,max m,mm ^ (5) 



m,max mi 



where subscript i represents pixel i. From Eq.(5) we can see that /? can be obtained conveniently from 
the data shown on Figure. 1 . 

2.4 Estimation of Other Core Variables 

After T s , T v , fi s and/? v are obtained using the above methods, net radiation at the soil surface and 
at the vegetation surface (i? n , soil, ^n,veg) can be calculated following Eqs. (6) and (7), and thereby soil 
heat flux and net radiation of mixed pixel can be estimated using the empirical formulation of Eq. (8) 
and Eq. (9), as used widely in previous studies [15-17]. 

R n,soil = ^0 (1 _ & soil ) + aS skyT s A ky ~ aS soil^toil ( 6 ) 
R n,veg ~ ^0 (1 — ®veg ) ^sky^sky ~ ^veg^veg (^) 

G » 0.3(1 - 0.9 f)R n soil (8) 

R n =fRn, S oU+(l-f)Rn, Veg (9) 

where So is the solar incident total radiation and is regarded as spatially uniform for clear sky 
conditions at the regional scale, usually obtained from standard meteorological station; a so ii and a ve g 
are the albedo of bare soil and vegetation, also calculated by the PCACA algorithm (seen from 
Appendix H); e sky is the average sky emissivity and is approximately set to 1.0 in the study. 7^ is 
average sky temperature and usually approximates to the temperature at 37° view sky angle [18], 
detailed calculation about which will be described in the Section 3. 

In terms of the Bowen-ratio energy balance method, soil evaporation (AEsoii) and vegetation 
transpiration (AEVeg) can be retrieved based on Eq. (10). When energy exchange between vegetation 
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and bare soil is neglected, XE of a mixed pixel can be described as a linear combination of XE so n and 
AEVeg, expressed as Eq. (11). 

^ = /^ veg +(i-/)^ o;7 (ii) 

3. Two improvements for the two-layer model 

J. i locating the true dry edge in T m -f space 

As mentioned above, the locations of true dry edge and the true wet edge are crucial in the 
application of PCACA algorithm. The previously used method to determine the trapezoid boundary 
often leads to uncertainties because of subjective judgement. To reduce the errors from this respect, a 
physically based method, which takes account of the surface energy balance, is presented in this study. 

According to Figure. 1, the four corner points, r s d, T sw , T Y & and JVwCan determine the envelop shape 
of the trapezoid, that is to say, as long as their values are obtained. Consequently, the true dry edge and 
the true wet edge can be determined. In the study, surface energy balance method was adopted to 
compute their values. For pixels at the true dry edge, Eq.(12) is used because XE=0. Substituting R n , G, 
H for Eq.(9), Eq.(8) and Eq.(13), we obtained Eq.(14) for T s d. In the same way, JVd is formulated as 
Eq.(15). 

R n -G = H (12) 
H = P Cp(T s -T a ) (13) 



r 

'a 



0.7[5 0 (1 - CCsd ) + ^skyTsky ] + — ^ 
T _ r sda (14) 

1 sd ~ 



'sda 



0.7[S 0 (1 -a vd ) + as sky TX ] + ^T vda (15) 



vda 



pCp 



r. 



vda 



3 

veg 1 vd 



where the subscripts sd, vd represent true driest bare soil and true driest full-cover vegetation, 
respectively, p is density of air, c p is the volumetric heat capacity of air, r a is the air dynamic 

resistance, a is the S-B coefficient, « S d is the albedo of dry bare soil. 

From these two equations, we can see that r s d or T v & could not be iteratively computed until 
parameters of a s d, a v d, T^ y , r s d a , r v d a , r s d a , 7Vd a are acquired. From the T m -f space, the observed driest 
bare soil and the observed driest full-cover vegetation can be found. Although there are surface 
temperature differences between them and the true driest bare soil and true driest full-cover vegetation, 
they represent relative driest and wettest soil conditions. It means that aerodynamic resistance and air 
temperature of the observed driest bare soil and the observed driest full-cover vegetation are 
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approximate to that of T s &, 7Vd points under the conditions that the variations in atmospheric conditions 
are very small due to the spatially uniformity. In the study, we chose 50 pixels around the upper-left 
corner in the trapezoid representing the observed driest bare soil and 50 pixels around upper-right 
corner in the trapezoid representing the observed driest full-cover vegetation, and the highest T a and r a 
in each 50 pixels were selected as T s d a and r s d a , TVda and r v d a , respectively. The method of retrieving the 
spatial distribution of T a and r a will be described in the following. The calculation of r a requires e a and 
r s , the retrievals about which will also be described. As for a sd and a v d, we selected the lowest albedo 
values of bare soil and full-cover vegetation within the whole scene of remotely sensed image (pixels) 
because low albedo would result in high surface temperature at the same soil water content, judging 
from Eq. (14) and Eq. (15). The item of cr£,T s t, represents downward long wave radiation that usually 
has small spatial variability for clear sky, thus can be obtained from the measurements of 
meteorological stations at the satellite overpass time. The following are the calculations of spatial 
distributions of T a , e a , r s and r a \ 



a) Estimation of air temperature (T a ) 



Surface temperature, as a heat or cold source, influences the variations of air temperature by heating 
or cooling near-surface atmosphere. In most cases, high surface temperature is accompanied by high 
air temperature and low surface temperature is accompanied by low air temperature. By using this 
relationship between them and assuming the following ratio expressed in Eq. (16), an interpolation 
method is presented to map the distribution of air temperature. The numerator and the denominator 
both represent radiant emission. Because the effects of neighborhood pixels diminish as the distance to 
the pixel increases, a weight (w) derived from an inverse distance weight (IDW) method is assigned to 
each neighboring observations, and is given in Eq. (17). 



n 

'ai ' 1 ai _ „ ^4 



i=\ _ " "aj ^aj 



°Sai< (16) 



T 4 

4 (7 • S m j ■ 1 m j 

mi - 1 mi 



(=1 



^='-24<^)" 7 <17) 

where i represents the pixel where air temperature measurement is made, j indicates the pixel where no 
measurements are available and estimate is required; n is the number of air temperature observations; d 
is the distance between pixel i and pixel j; p is an exponent. The higher the exponent is, the larger the 
influence of the closest observations on estimate is; p is set to 2 in the study; a is the Stephen 
Boltzmann constant; e a is the emissivity of air calculated from Eq. (17) [19]. 



b) Estimation of actual vapor pressure near surface (e a ). 



Like the relationship between air temperature and surface temperature, there are also strong 
interactions between near-surface actual vapor pressure and soil moisture. If there is no horizontal and 
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vertical advection, vapor in the atmosphere mainly comes from soil water by an evapotranspiration 
process. Contrarily, the vapor gradient between surface and air influences the intensity of 
evapotranspiration. By assuming the following ratio relationship between near-surface actual vapor 
pressure and soil water, we interpolated e a , where soil moisture status was characterized quantitatively 
by the soil apparent thermal inertia: 

Rn 

2(T-T. )Af 2 

V $ nun / 

n 

H W i' e ai 

i=l e aj (18) 

<f R n 

2^ w i 



nj 



,=] (T mi T m i n ) (T m j T mm ) 



where i, j, W\, T s have the same meanings as in Eqs. (16) and (17); T min is the minimum surface 
temperature during the daytime which usually occurs before sunrise when R n =0 and for all pixels it can 
be assumed to be the same value [20]. At is the time interval between sunrise and the satellite overpass 
time. 



c) Estimation of surface resistance to evapotranspiration (r s ) 



In the retrieval of evapotranspiration, r s is used to correct the difference between the vapor pressure 
at the surface (e s ) and the saturated vapor pressure at the evaporating front (e s *). In theory, r s ranges 
from 00 to 0 corresponding to the surface conditions of potential evapotranspiration and zero 
evapotranspiration, respectively, namely the true dry edge and the true wet edge. However, in reality 
the condition of zero evapotranspiration (r smax = °°) rarely occurs for vegetated surface even in semi- 
arid environment primary due to root zone soil water uptake, consequently, we selected a pixel closest 
to the observed driest bare soil where a meteorological station is located to calculate r smax by IE, e a , u 
measurements at the satellite overpass time. r smax is about 1000 (s m" 1 ) according to the calculation, 
which is in agreement with Qiu's observations [21]; r sm i n is set to 0 in the study. 

After the upper and lower limits of r s are determined, r s is interpolated linearly within each / 
interval between the lowest and the highest temperature. Following the illustrations in Figure 1, r s ; for 
a pixel at (fi,T\) equals: 

T si- T siMn , . ( 19 ) 

si rp rp V'.smax 's mm/ 1 'smm ' 

Asz',max — * si,mm 

The interpretation is similar to the method used by Stisen et al. [11]. The difference is that Stisen et al. 
used Priestly-Taylor parameter to represent an effective surface resistance to evapotranspiration. 



d) Estimation of aerodynamic resistance (r a ) 



Besides air temperature, aerodynamic resistance is a site-specific variable and can not be retrieved 
directly by remote sensing. Although Monin-Obukhov similarity theory has been widely used to 
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estimate it, the accurate calculations for the spatial distributions of roughness length (Z 0 ), wind speed 
(u) and the atmospheric stability parameters are very difficult. In this study, we adopted the energy 
balance method. 

An expression of r a is obtained on the basis of the energy-balance by substituting Rn, G, H and XE 
for Eq. (9), Eq. (8), Eq. (13) and Eq. (20): 



* 



+ a 2 r a + a 3 = 0 



XE= P C p( e s e a) (20) 
r( r s+ r a) 

ai =(0J + 0.27f)- r -R n (21) 
a 2 = (0.7 + 0.27 f)/ -R„-r s - pc p y{T m -T a )- pc p {e s - e a ) 
a 3 =-pc p Y(T m -T a )-r s 

where y is the psychometric constant, e s * is the saturated vapor pressure at T m estimated using the 
classic formulation with regard to surface temperature, as Eq. (22), other parameters have the same 
meanings as the above illustrations. Here, T a , e a and r s are estimated by the approaches mentioned 
previously. 

e;=6.1078exp(i 727 ^) (22) 
T m + 237.3 

Above all, all necessary parameters to calculate r s d, 7Vd can be retrieved according to the Eqs. (12) 
- (22) in combination with limited ancillary data. Here, T m and a m is obtained from MODIS (Moderate 
Resolution Imaging Spectroradiometer) standard land data products, MOD 11 and MOD02 [22, 23], 
through the NASA Earth Observing System Data Gateway. 



3.2 Locating the true wet edge in T m -f space 



Many studies [9, 24, 25] indicated that the dry edge is more evident than the wet edge in the T m -f 
space. There often exist outlying points exhibiting a tail toward low values of temperature and / Such 
points may represent anomalous surfaces mainly due to cloud contamination and usually are discarded 
from the analysis [10]. Thus, although r sw and ^also can be parameterized on the basis of energy 
balance equation like and 7Vd, the inputs (a s d, a v d, ^sda, ^vda, T'sda, TVda) cannot be obtained by the 
above methods. Taking account of the fact that the surface radiant temperature of dense vegetation is 
very close to the ambient air temperature [9, 26, 27], we took the average air temperature at 7=1 as r vw . 
As for T sw , we adopted an approximation of using the surface temperature of standing water body 
(such lake) as the surface temperature of true wet bare soil J" sw , that is to say, standing water body is 
regarded as the surface of potential evapotranspiration. In fact, it is not uncommon to find some 
patches of standing water body in a remote sensed image. In this study, the surface temperature of 
Dongping Lake (35.965°, 116.81°; water area: 209 km 2 ) was used as T sw . In applications, we found 
that the pixels with mixed surface temperatures below the true wet edge all scattered around the cloud 
pixels and the coastline pixels. 

3.3 Physical illustrations for the uncertainties using the above locating methods 
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Using the above methods, the positions of the true dry edge and the true wet edge in the triangle can 
be located. When surface radiant temperature is mainly dominated by surface soil water content, the 
following relationship between surface temperature and soil water content is tenable (Figure 2), which 
has been illustrated in the Qiu's experiment [21]. Figure 2 suggests that when the soil water content 
ranges from 80% to 100% of the soil saturation, evapotranspiration will happen approximately at the 
level of potential Evapotranspiration. On the contrary, when the range of the soil water content is from 
the wilting point to the driest, evapotranspiration is close to zero, and in the two extreme conditions, 
surface temperature can be regarded as a constant due to the same evaporative cooling effect. On the 
basis of this, although the true driest and the true wettest points can't be determined using the above 
methods, the errors will be acceptable in the determination of T^, T sw , T v d and T vw due to the constant 
surface temperature. 



Figure 2. Relationship between surface temperature and soil water content. 




3.4 Elimination of the effects of atmospheric conditions on surface evapotranspiration 

Same as the triangle method, an important assumption for the PCACA algorithm is that the surface 
evapotranspiration is primarily constrained by soil water availability, based on which T so u values for 
all pixels with an equal water availability are identical, so are for T ye% values, that is to say, the 
configuration of T s and / is mainly caused by the variation of the soil water availability and is 
irrelevant to the differences in atmospheric conditions and surface properties [11]. However, in reality 
it is not true that the atmospheric conditions and surface properties are entirely homogeneous within 
the large study area. There exist some differences eg. in wind speed, air temperature, surface 
roughness length and surface albedo, which are all site-specific variables. As a result, it will lead to the 
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the variability of the surface evapotranspiration and thereby influence the surface temperature and the 
configuration of T s and / To satisfy the assumption, the effects of these factors need to be eliminated. 

In the study, the effects of four controlling factors (r a , e a , T a and albedo) on surface temperature 
were eliminated. According to the assumption, r a , e a , T a and albedo values of each pixel should be 
equal, namely they are spatially homogeneous, to ensure that only water availability controls surface 
temperature. To meet this requirement, the average values of r a , e a , T a and albedo in the image are 
assigned to each pixel and a new mixed surface temperature (T m i) is calculated based on energy 
balance equation. T m i is the assumed temperature controlled only by soil water availability and 
therefore the new trapezoid constructed by T m and / is more meaningful for the PCACA algorithm. It 
has to be noted that although the configuration of T m -f space is different from that of the T m -f space, 
the locations of the true dry edge and the true wet edge in the trapezoid don't change since the two 
theoretical extreme soil water status at the satellite overpass time are the same. On the basis of T m , 
Tveg (vegetation temperature controlled only by soil water availability) and T so n (soil temperature 
controlled only by soil water availability) can be obtained. Now the problem is that the calculated T veg 
and Tsoii cannot represent the actual temperature of vegetation and soil after the above transformation., 
Therefore they can't be directly used in the consequent calculations. To estimate the true vegetation 
temperature and true soil temperature (T veg and T so n) from T veg and T so n , the following method is used. 

Assuming that the thermal energy fraction assigned to the vegetation and the soil are constants, the 
following expression [Eq. (23)] can be approximated. The left item represents the proportion of 
thermal energy difference for vegetation or soil to the total thermal energy difference, and the right 
item represents the proportion of vegetative or soil thermal energy to the total thermal energy. It 
suggests that changes of environmental conditions bring same effects on vegetation temperature, soil 
temperature and mixed surface temperature, for example, T so n, T veg and T m values would all increase 
with the increase of solar radiation intensity and all influenced by wind speed. It is not likely to happen 
that the T^u increases, but T veg decreases in the same atmospheric conditions. Therefore, the 
assumption is reasonable in practice and it is different from Beer's law using fractional vegetation 
cover as the weight to partition thermal energy. Using this equation, T so n 4 - T veg can be solved as Eq. 
(24). 

T 4 —T' 4 T 4 t 4 —T' 4 T 4 

veg 1 veg ^ 1 veg 1 soil 1 soil „ 1 soil (23) 

rp4 _ rp'4 ~ rp4 ' J>4 _ rp '4 ~ rp 4 ^ ' 

mm m mm m 



T 4 

T 4 -T 4 =-S-(T' 4 -T' 4 (24) 

veg .soil rp '4 V veg soil ) V ' 

m 

Combined with Eqs. (1) and (24), 7" so ii and r veg can be solved. r so ;i and r veg can be directly obtained 
from Eq. (23), at the same time the uncertainties of the solution can be reduced to the minimum by this 
method because the difference between T so n and T veg is used instead of the absolute temperatures and 
Eq.(l) helps to maintain the thermal energy balance in the calculation of T so n and T yes values. 
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4. Study area and field measurements 

The study area is located in the North China Plain and ranges from 35.2N to 40.84N in latitude, 
from 113.68E to 119.54E in longitude. The land use in the area is dominated by the rotating cropping 
of winter wheat and summer maize. Millet, soybean and cotton are also scattered planted in summer 
[28]. According to the traditional tillage practice, winter wheat is sown in early October, harvested in 
early or mid June next year, and summer maize is planted in early to mid June and harvested at the end 
of September. The soil is mostly silt, light loam and medium loam. Annual precipitation is about 600 
mm, more than 50% of which falls during the summer monsoon between July and September. The 
groundwater table varies from 1.5 m to 3.5 m with an average of 2.5 m. 

In the study, field measurements from 135 standard meteorological stations were used. The 
measurements include air temperature and actual vapor pressure at 2 m height above the surface, solar 
incoming radiation, surface radiative temperature, wind speed, upward longwave radiation, upward 
shortwave radiation, downward longwave radiation and downward shortwave radiation. Figure. 3 
shows the spatial distribution of these stations in the study area, where the triangle symbol (A) marks 
the location of Yucheng Agro-ecosystem Station that can represents the largest agricultural area in the 
North China Plain [29]. Sensible heat flux and latent heat flux have been continuously measured by the 
Eddy Correlation (EC) system which is composed of a 3D sonic anemometer and an open path 
CO2/H2O analyzer since 2002. H, XE and CO2 fluxes were originally sampled at 10 Hz and values 
averaged over 30 min were used in the study to validate the above mentioned methods. The rectangle 
symbol (■) shows the location of Dongping Lake. Water surface temperature has been measured for 
seven years since 2001 and was used to validate the MODIS land surface temperature products used in 
this study, and thereby to ensure the accuracy of r sw values. 

Figure 3. Spatial distribution of standard meteorological stations in the study area. 
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5. Satellite data 

MODIS land data products, including MOD11 (Land Surface Temperature), MOD03 (Geolocation 
Data Set), MOD05 (Total Precipitable Water), MOD02 (Calibrated Geolocated Radiance) and MOD35 
(Cloud Mask), for clear sky during springtime between Mar and June when winter wheat is the 
dominant crop were used in the study. The visible channels of MOD02 were processed with 
atmospheric correction using the SMAC algorithm [30] with the combination of the MOD03 and 
MOD05 data. The cloud was masked out based on MOD35 data. The basic variables used as inputs in 
the model consist of albedo, fractional vegetation cover and surface temperature, based on which the 
calculation of R n , PCACA algorithm and Layered Energy-separating Algorithm were applied. 

The instantaneous surface albedo was obtained by averaging reflectance values for several visible 
and near-infrared channels with the wavelength as the weight. This would introduce some errors 
because channel reflectance values observed by a satellite are reflectance at only one Sun-target-sensor 
configurations and it is generally different from the hemispherical reflectance. But we assumed its 
influence to be small as pointed by Nishida [27]. The fractional vegetation cover was estimated from 
the normalized difference vegetation index (NDVT) [27], as Eq.(25). 

NDVI- NDVI mtn 



f = 



NDVI max -NDVI mm 



(25) 



where NDVI max and NDVI m m are NDVI values for full cover vegetation (f=\) and bare soil (f=0). 
According to the field NDVI measurements performed at Yucheng station, NDVI m i n =0.09 and 
NDVI max =0JS. Other variables, such as R n , G, T a , e a , r s , were obtained using the above mentioned 
methods. 

6. Results and Discussion 

Figure 4 shows the comparisons between T soi \, T veg obtained from the original model and T soi \ , r veg 
obtained from the improved model at an iso-line of equal water availability. 

Figure 4. Comparisons between T so n, T veg obtained from the original model and T so n , r veg 
obtained from the improved model at iso-line of equal water availability. 
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Apparently, the results of T so n , T veg are in closer agreement with the above assumptions that soil 
surface temperature for all pixels at an iso-line are identical, so are the vegetative surface temperature. 
It provides a better physical foundation for the model. 

Due to the absence separate measurement of the soil evaporation, soil heat flux, vegetation 
transpiration and vegetation heat flux from the field sit, estimates of R n -G and XE from the model were 
compared to the measurements from Yucheng station, as shown in Figure 5 for available energy (R n - 
G), in Figure 6 for XE. 

Figure 5. Modeled versus measured available energy, R Q -G. The line represents a 1: 1 
relationship. 




Figure 6. Modeled versus measured latent heat flux, XE. The line represents a 1: 1 
relationship. 
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Because direct comparison of /LE so ii and XE veg cannot be performed, the relationship between 
modeled /LE veg and measured CO2 were used to indirectly validate the two-layer model since CO2 
fluxes are closely related to /LEVeg [31, 32]. Figure 7 shows the scatter plot between modeled /LEVeg and 
measured CO2 flux. Note that the minus sign (-) for CO2 fluxes means that the flux transferring 
direction is from up to down. The performance of the model was evaluated using the root mean 
squares difference (RMSD) and the mean absolute difference (MAD), which are defined as Eq.(26) 
and Eq.(27), respectively. 

n 

RMSD = [^(P i -O i ) 2 /n] 1 ' 2 (26) 
1=1 

n 

MAD = Y\Pi-Oi\ln ( 27 ) 

z=l 

where n is the number of observations, P represents the model estimated value, O represents the 
observed value. 

Figure 7. Seasonal variation of the modeled /LE veg and measured CO2 flux during winter 
wheat growing period 2005 and 2006 (from March to June). 
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The figures demonstrate that (1) no significant bias is found between the modeled and measured R n - 
G with RMSD=46.3 and MAD=40.2. The high correlation coefficient (r 2 =0.79) suggests that soil 
and i? n , ve g calculated from T so n and T veg are reasonable because Rn is calculated by them; (2) estimates 
of XE tends to overestimate for lower XE values and shows larger bias than R n -G resulting in 54.1 of 
RMSD and 47.7 of MAD, which implies a larger uncertainty of the model in low XE conditions; (3) in 
general, seasonal variations in modeled XE veg and measured C0 2 flux shows good agreement though 
few points in 2006 showed different variations perhaps influenced by horizontal/vertical advection or 
other factors. Vegetation transpiration generally increases with the crop growth, at the same time, since 
vegetation absorbs more CO2 due to the more active photosynthesis in the daytime, CO2 fluxes also 
increases. On the contrary, XE weg and CO2 fluxes both became very small after crop harvest or in winter. 
In the study, winter wheat was harvest in June, therefore XE veg and CO2 fluxes rapidly decreased on 
June 15 and June 19 in 2005. Figure 7 indirectly proves the soundness of the two-layer model to 
estimate the soil evaporation and vegetation transpiration separately. 

In fact, there are two other factors which contribute to the uncertainty in the above validation. One 
is that point measurements usually can not represent the whole MODIS pixel (1 km 2 ) because of the 
large scale disparity between them. The other is that the observed R n , XE, H and G values at field 
cannot meet the surface energy balance closure resulting from many possibilities, such as instrumental 
errors, horizontal/vertical advections [33-35], however, the basis of the model used in the study is the 
surface energy balance, as most remote sensing models of evapotranspiration, such as SEBS [36], 
SEBAL [37], N95 [3]. Therefore, it is not uncommon to see the differences between modeled fluxes 
and observed ones. Compared with other studies, the RMSD of 54.1 (w/m 2 ) for modeled^ are in a 
acceptable range [17, 27, 38, 39]. 

Figure 8. XE veg maps on May 2, 2005 in North China Plain retrieved by the improved 
model (left image) and the original model (right image). 
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Figure 8 shows the /LEVeg maps on May 2, 2005 for the North China Plain retrieved by the improved 
model and the original model, respectively. Differences were found between them. The /LEVeg value 
from the improved model varies from 222 (w/m 2 ) to 456 (w/m 2 ), while the original /LEVeg values varied 

2 2 

from 50 (w/m ) to 560 (w/m ). In terms of the experiences and the knowledge on the vegetation 
transpiration, a small /LEVeg dynamical range is more reasonable in the rapid growing season of winter 
wheat, that is to say, the results from the improved model is better than that of the original model 
although it is difficult to quantitatively evaluate the two results due to the absence of field 
measurements of the vegetation transpiration and soil evaporation at the satellite pixel scale. 
Furthermore, according to the above illustrations, the physical basis of the improved model is 
strengthened against the original model. 

7. Conclusions 

This paper presented two improvements of a two-layer model for the estimation of the land surface 
heat fluxes. The weakness in the original model was identified: (1) a subjective method to determine 
the true boundary lines from the scatter plot for the surface temperature of mixed pixel versus the 
fractional vegetation cover; (2) the assumption that the configuration of T m — f space is mainly 
controlled by soil water availability, are not physically realistic. To investigate the two issues and 
obtain more realistic separation of the surface temperature for the soil and vegetation components, and 
thereby more accurate partitioning of the surface fluxes for the soil and vegetation components, 
introducing the surface energy balance method solves the problem (1) and the method of eliminating 
the effects of four controlling factors (r a , e a , T a and albedo) on surface temperature solves the problem 
(2). At the same time, the interpolation methods of T a and e a , the methods for acquiring spatial 
distributions of r a and r s were also investigated. Finally, the improved model was applied to the North 
China Plain. The results showed good agreement with in situ terrestrial surface fluxes measurements. 
Furthermore, by comparing the seasonal variations of vegetation transpiration and C0 2 flux, and the 
vegetation transpirations retrieved by the improved model and the original model, the effectiveness of 
the improved two-layer model for estimating soil evaporation and vegetation transpiration were 
indirectly proved. 

Acknowledgements 

This work was supported jointly by the State Key Development Program for Basic Research of 
China with grant number 2007CB714401-3, the Knowledge Innovation Project of IGSNRR, CAS 
(grant No. GXIOG-A05-1 1), the Knowledge Innovation Program of CAS (grant No. KZCX2-YW- 
326) , the Young Talents Forefront Project of IGSNRR, CAS (grant No.07v70050SZ), the National 
Key Project of Scientific and Technical Supporting Programs Funded by Ministry of Science & 
Technology of China (No. 2006BAC08B0407) and the Program of "One Hundred Talented People" of 
the Chinese Academy of Sciences (CAS) 



Sensors 2008, 8 



6183 



Appendix I 



The relationship between the temperature of mixed pixel and vegetative fractional cover for 
8~14//m can be expressed as: 

VZrJm = ^vegfTveg + ° £ soil( l ~ f) T soil ( AU ) 

where T m , T veg and T so n are respectively the surface temperature of mixed pixel, vegetation canopy and 
soil surface. e m , e veg and £ so ii are the emissivities of mixed pixel, vegetation and bare soil; c is the 
Stefan-Boltzmann constant. / is fractional vegetation cover of mixed pixel. Evidently, r veg and r so ii 
cannot be retrieved only using this equation. Computing the differential coefficient of T m to / and 
assuming 

df df 

we can get Eq. (AI.2): 

4 r 3^k + T 4d£ ]IL= T 4 _ t a (AI.2) 

"mm m ^ veg veg "soil soil 

Multiplying/in the two sides of Eq (AI.2) and rearranging this equation, we obtained: 
f(A P T 3 ^ m i T 4 ^ Sm \ - fr T A - fp T A CAI 3) 

J \ to m ± m m ^ ) J veg veg J soil soil v " ' 

Integrating Eq. (AI.2) with Eq. (AI.3), T^u can be expressed as: 

Tson = —Vsjl - f(T: ^ + *sjl ^)] (AI.4) 



in the same way, 7V eg can be expressed as: 

T ves = \ —leJm + (1 - Wm + ^Jl ^)] \ ( AL5 > 



Sveg df df 



1/4 



ds ds 
= m the study, 



According to Eq. (3), = s veg - s soil . In the study, = 0.02 due to the fixed values of e V e g =0.97 



ds 

and e so ii=0.95. In terms of the numerical simulations, the small — — value has small influence on the 

df 

difference between T soi \ and T veg , therefore we ignored the difference between e veg and e SO ii. This 
simplification means that Eq. (AI.l) was simplified as Eq. (AI.6), which is also reasonable in 
applications [5]. 

t=^;+a-/)C (ai.6) 

Furthermore, using Eq.( AI.5) minus Eq.( AI.4), the following expression can be simplified as: 

df 



T 4 -T 4 , *4T 3 ^- (AI.7) 

veg soil m ■, s- 
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In fact, T* - Tsoii ~ ^(T veg - T soil )T^ when no large temperature difference between vegetation surface 



dT 

and soil surface occurs, so — — can be formulated as: 

df 



dTm-r _ T , (AI.8) 

^j. veg ^soil ' 



Appendix H 



Albedo of mixed pixel is the weighed-average value of soil albedo and vegetation albedo, showed 

as Eq. (AH.l). Computing the differential coefficient of a m to / and assuming ^ a veg _ q da soil _ q ? we 

df ' df 

obtained Eq.( AH.2). 

a m =(l-f)a soil +fa veg (AH.l) 
^ = a -a , (AH.2) 

Integrating Eq. (AH.l) with Eq. (AH.2), a so ii and a veg can be expressed as: 

*son=* m -f d -^ ( AH - 3 > 
df 

a veg =<x m +(l-f)^ (AH.4) 

Evidently, as the method of decomposing mixed surface temperature T m , PCACA also can be applied 
to separate albedo of mixed pixel using the configuration of a m - f space. 
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